#packages
library(haven)
library(lattice)
library(misty)
library(Matrix)
library(lme4)
library(msm)
library(sandwich)
library(sjPlot)
library(lmerTest)
library(backports)     
library(carData)
library(effects)       
library(ggplot2)       
library(interactions)  
library(psych)         
library(plyr)
library(ggpubr)


#################################################################################################################
#Prepare data for covariates models
########################################################################################################################

dataset_rmssdco <- dataset[ , c("id", "resting_rmssd", "reactivity_rmssd", "resting_met", "reactivity_met", "resting_supine", "reactivity_supine", "mvpa", "bmi", "chronicstress", "conttime", "caffeine_full", "nicotine_full", "alcohol_full", "ambientnoise_full", "anticipatedstress_full")]
dataset_rmssdco <- dataset_rmssdco[complete.cases(dataset_rmssdco),]

#RMSSD reactivity 
dataset_rmssdco$rmssd_reactivity <- (dataset_rmssdco$reactivity_rmssd - dataset_rmssdco$resting_rmssd)

#Log transform RMSSD parameters
#RMSSD reactivity
dataset_rmssdco$reactivity_rmssd_log <- log(dataset_rmssdco$reactivity_rmssd)
#RMSSD resting
dataset_rmssdco$resting_rmssd_log <- log(dataset_rmssdco$resting_rmssd)

#supine position in minutes
dataset_rmssdco$resting_supine_minutes <- (dataset_rmssdco$resting_supine)*5
dataset_rmssdco$reactivity_supine_minutes <- (dataset_rmssdco$reactivity_supine)*5

#Scale independent variables
#Resting MET
dataset_rmssdco$resting_met_cent <- scale(dataset_rmssdco$resting_met, scale = FALSE)
#MVPA
dataset_rmssdco$mvpa_cent <- scale(dataset_rmssdco$mvpa, scale = FALSE)
#BMI
dataset_rmssdco$bmi_cent <- scale(dataset_rmssdco$bmi, scale = FALSE)
#Chronicstress
dataset_rmssdco$chronicstress_cent <- scale(dataset_rmssdco$chronicstress, scale = TRUE)
#RMSSD Resting
dataset_rmssdco$resting_rmssd_log_cent <- center(dataset_rmssdco$resting_rmssd_log, type = c("CWC"), cluster = dataset_rmssdco$id, value = NULL, as.na = NULL,
                                                   check = TRUE)
#Reactivity MET
dataset_rmssdco$reactivity_met_cent <- scale(dataset_rmssdco$reactivity_met, scale = FALSE)
#Conttime_hour
dataset_rmssdco$conttime_cent <- scale(dataset_rmssdco$conttime, scale = FALSE)

#outliers
#Reactivity RMSSD
dataset_rmssdco$rmssd_reactivity[which(dataset_rmssdco$rmssd_reactivity < (mean(dataset_rmssdco$rmssd_reactivity) - 3*sd(dataset_rmssdco$rmssd_reactivity)))] <- (mean(dataset_rmssdco$rmssd_reactivity) - 3*sd(dataset_rmssdco$rmssd_reactivity))
dataset_rmssdco$rmssd_reactivity[which(dataset_rmssdco$rmssd_reactivity > (mean(dataset_rmssdco$rmssd_reactivity) + 3*sd(dataset_rmssdco$rmssd_reactivity)))] <- (mean(dataset_rmssdco$rmssd_reactivity) + 3*sd(dataset_rmssdco$rmssd_reactivity))
#RMSSD resting
dataset_rmssdco$resting_rmssd_log_cent[which(dataset_rmssdco$resting_rmssd_log_cent < (mean(dataset_rmssdco$resting_rmssd_log_cent) - 3*sd(dataset_rmssdco$resting_rmssd_log_cent)))] <- (mean(dataset_rmssdco$resting_rmssd_log_cent) - 3*sd(dataset_rmssdco$resting_rmssd_log_cent))
dataset_rmssdco$resting_rmssd_log_cent[which(dataset_rmssdco$resting_rmssd_log_cent > (mean(dataset_rmssdco$resting_rmssd_log_cent) + 3*sd(dataset_rmssdco$resting_rmssd_log_cent)))] <- (mean(dataset_rmssdco$resting_rmssd_log_cent) + 3*sd(dataset_rmssdco$resting_rmssd_log_cent))
#RMSSD reactivity
dataset_rmssdco$reactivity_rmssd_log[which(dataset_rmssdco$reactivity_rmssd_log < (mean(dataset_rmssdco$reactivity_rmssd_log) - 3*sd(dataset_rmssdco$reactivity_rmssd_log)))] <- (mean(dataset_rmssdco$reactivity_rmssd_log) - 3*sd(dataset_rmssdco$reactivity_rmssd_log))
dataset_rmssdco$reactivity_rmssd_log[which(dataset_rmssdco$reactivity_rmssd_log > (mean(dataset_rmssdco$reactivity_rmssd_log) + 3*sd(dataset_rmssdco$reactivity_rmssd_log)))] <- (mean(dataset_rmssdco$reactivity_rmssd_log) + 3*sd(dataset_rmssdco$reactivity_rmssd_log))
#Reactivity MET
dataset_rmssdco$reactivity_met_cent[which(dataset_rmssdco$reactivity_met_cent < (mean(dataset_rmssdco$reactivity_met_cent) - 3*sd(dataset_rmssdco$reactivity_met_cent)))] <- (mean(dataset_rmssdco$reactivity_met_cent) - 3*sd(dataset_rmssdco$reactivity_met_cent))
dataset_rmssdco$reactivity_met_cent[which(dataset_rmssdco$reactivity_met_cent > (mean(dataset_rmssdco$reactivity_met_cent) + 3*sd(dataset_rmssdco$reactivity_met_cent)))] <- (mean(dataset_rmssdco$reactivity_met_cent) + 3*sd(dataset_rmssdco$reactivity_met_cent))
#Resting MET
dataset_rmssdco$resting_met_cent[which(dataset_rmssdco$resting_met_cent < (mean(dataset_rmssdco$resting_met_cent) - 3*sd(dataset_rmssdco$resting_met_cent)))] <- (mean(dataset_rmssdco$resting_met_cent) - 3*sd(dataset_rmssdco$resting_met_cent))
dataset_rmssdco$resting_met_cent[which(dataset_rmssdco$resting_met_cent > (mean(dataset_rmssdco$resting_met_cent) + 3*sd(dataset_rmssdco$resting_met_cent)))] <- (mean(dataset_rmssdco$resting_met_cent) + 3*sd(dataset_rmssdco$resting_met_cent))
#Resting supine
dataset_rmssdco$resting_supine_minutes[which(dataset_rmssdco$resting_supine_minutes < (mean(dataset_rmssdco$resting_supine_minutes) - 3*sd(dataset_rmssdco$resting_supine_minutes)))] <- (mean(dataset_rmssdco$resting_supine_minutes) - 3*sd(dataset_rmssdco$resting_supine_minutes))
dataset_rmssdco$resting_supine_minutes[which(dataset_rmssdco$resting_supine_minutes > (mean(dataset_rmssdco$resting_supine_minutes) + 3*sd(dataset_rmssdco$resting_supine_minutes)))] <- (mean(dataset_rmssdco$resting_supine_minutes) + 3*sd(dataset_rmssdco$resting_supine_minutes))
#Reactivity supine
dataset_rmssdco$reactivity_supine_minutes[which(dataset_rmssdco$reactivity_supine_minutes < (mean(dataset_rmssdco$reactivity_supine_minutes) - 3*sd(dataset_rmssdco$reactivity_supine_minutes)))] <- (mean(dataset_rmssdco$reactivity_supine_minutes) - 3*sd(dataset_rmssdco$reactivity_supine_minutes))
dataset_rmssdco$reactivity_supine_minutes[which(dataset_rmssdco$reactivity_supine_minutes > (mean(dataset_rmssdco$reactivity_supine_minutes) + 3*sd(dataset_rmssdco$reactivity_supine_minutes)))] <- (mean(dataset_rmssdco$reactivity_supine_minutes) + 3*sd(dataset_rmssdco$reactivity_supine_minutes))
#MVPA
dataset_rmssdco$mvpa_cent[which(dataset_rmssdco$mvpa_cent < (mean(dataset_rmssdco$mvpa_cent) - 3*sd(dataset_rmssdco$mvpa_cent)))] <- (mean(dataset_rmssdco$mvpa_cent) - 3*sd(dataset_rmssdco$mvpa_cent))
dataset_rmssdco$mvpa_cent[which(dataset_rmssdco$mvpa_cent > (mean(dataset_rmssdco$mvpa_cent) + 3*sd(dataset_rmssdco$mvpa_cent)))] <- (mean(dataset_rmssdco$mvpa_cent) + 3*sd(dataset_rmssdco$mvpa_cent))
#BMI
dataset_rmssdco$bmi_cent[which(dataset_rmssdco$bmi_cent < (mean(dataset_rmssdco$bmi_cent) - 3*sd(dataset_rmssdco$bmi_cent)))] <- (mean(dataset_rmssdco$bmi_cent) - 3*sd(dataset_rmssdco$bmi_cent))
dataset_rmssdco$bmi_cent[which(dataset_rmssdco$bmi_cent > (mean(dataset_rmssdco$bmi_cent) + 3*sd(dataset_rmssdco$bmi_cent)))] <- (mean(dataset_rmssdco$bmi_cent) + 3*sd(dataset_rmssdco$bmi_cent))
#Chronic stress
dataset_rmssdco$chronicstress_cent[which(dataset_rmssdco$chronicstress_cent < (mean(dataset_rmssdco$chronicstress_cent) - 3*sd(dataset_rmssdco$chronicstress_cent)))] <- (mean(dataset_rmssdco$chronicstress_cent) - 3*sd(dataset_rmssdco$chronicstress_cent))
dataset_rmssdco$chronicstress_cent[which(dataset_rmssdco$chronicstress_cent > (mean(dataset_rmssdco$chronicstress_cent) + 3*sd(dataset_rmssdco$chronicstress_cent)))] <- (mean(dataset_rmssdco$chronicstress_cent) + 3*sd(dataset_rmssdco$chronicstress_cent))

#Rescale
dataset_rmssdco$conttime_cent <- dataset_rmssdco$conttime_cent/24
dataset_rmssdco$mvpa_cent <- dataset_rmssdco$mvpa_cent/60

################################################################################################################
#Model with covariates
########################################################################################################################

model_covariates <- lmerTest::lmer(rmssd_reactivity ~ resting_rmssd_log_cent + resting_met_cent + reactivity_met_cent + resting_rmssd_log_cent*reactivity_met_cent + resting_supine_minutes + reactivity_supine_minutes + conttime_cent + caffeine_full + alcohol_full + anticipatedstress_full + ambientnoise_full + mvpa_cent + bmi_cent + chronicstress_cent + (1|id),data = dataset_rmssdco)
tab_model(model_covariates, show.se = TRUE)

################################################################################################################
#Figures
########################################################################################################################

dataset_rmssdco$predicted <- predict(model_covariates)

subsetlowmet <- subset(dataset_rmssdco, reactivity_met<3)
subsetmediumhilfe <- subset(dataset_rmssdco, reactivity_met>3)
subsetmediummet <- subset(subsetmediumhilfe, reactivity_met<6)
subsethigmet <- subset(dataset_rmssdco, reactivity_met>6)

plotlight <- ggplot(data=subsetlowmet, aes(x=resting_rmssd_log_cent, y=predicted, group=factor(id), colour="gray"), legend=FALSE) +
  geom_smooth(method=lm, se=FALSE, fullrange=TRUE, lty=1, size=.5, color="gray40") +
  geom_smooth(aes(group=1), method=lm, se=FALSE, fullrange=FALSE, lty=1, size=2, color="blue") +
  xlab("Resting lnRMSSD") + ylab("Predicted RMSSD Reactivity (ms)") +
  theme_classic() +
  theme(axis.title=element_text(size=18),
        axis.text=element_text(size=14),
        plot.title=element_text(size=18, hjust=.5))
plotlight <- plotlight + coord_cartesian(xlim = c(-1.5, 1.5), ylim = c(-40, 20))

plotmedium <- ggplot(data=subsetmediummet, aes(x=resting_rmssd_log_cent, y=predicted, group=factor(id), colour="gray"), legend=FALSE) +
  geom_smooth(method=lm, se=FALSE, fullrange=TRUE, lty=1, size=.5, color="gray40") +
  geom_smooth(aes(group=1), method=lm, se=FALSE, fullrange=FALSE, lty=1, size=2, color="orange") +
  xlab("Resting lnRMSSD") + ylab("Predicted RMSSD Reactivity (ms)") +
  theme_classic() +
  theme(axis.title=element_text(size=18),
        axis.text=element_text(size=14),
        plot.title=element_text(size=18, hjust=.5))
plotmedium <- plotmedium + coord_cartesian(xlim = c(-1.5, 1.5), ylim = c(-40, 20))

plothigh <- ggplot(data=subsethigmet, aes(x=resting_rmssd_log_cent, y=predicted, group=factor(id), colour="gray"), legend=FALSE) +
  geom_smooth(method=lm, se=FALSE, fullrange=TRUE, lty=1, size=.5, color="gray40") +
  geom_smooth(aes(group=1), method=lm, se=FALSE, fullrange=FALSE, lty=1, size=2, color="red") +
  xlab("Resting lnRMSSD") + ylab("Predicted RMSSD Reactivity (ms)") +
  theme_classic() +
  theme(axis.title=element_text(size=18),
        axis.text=element_text(size=14),
        plot.title=element_text(size=18, hjust=.5))
plothigh <- plothigh + coord_cartesian(xlim = c(-1.5, 1.5), ylim = c(-40, 20))

figure <- ggarrange(plotlight, plotmedium, plothigh,
                    labels = c("(A)", "(B)", "(C)"),
                    ncol = 3, nrow = 1)
show(figure)








